function [Z,V,etime,fval] = mcmc(X, Z, V, alphav, sigmax, sigmaz, sigmav, poisstrunc, algtype, loopi)
[N,K] = size(Z);
[D,tK] = size(V);
if K ~= tK
    error('Dimensions of Z & V must agree.');
end

if ~any(algtype == 1:4)
    error('cfg:invalidparam', ...
        ['Unknown ''algtype'': %d\n', ...
         '    1 for ''TYPE-I Gibbs sampling'';\n', ...          % TYPE-I: remove features in V whenever
         '    2 for ''TYPE-I MAP inference'';\n', ...           % their corresponding ones in Z get nullified
         '    3 for ''TYPE-II Gibbs sampling'';\n', ...         % TYPE-II: retain features in V even though
         '    4 for ''TYPE-II MAP inference''.'], algtype);     % their corresponding ones in Z get nullified
end

fval = [];
etime = 0;
wv = 2.^(N-1:-1:0);
for e = 1:1
    % sampling Z
    tStart = tic;
    newks = zeros(1,poisstrunc+1);
    m = sum(Z~=0);
    h = wv*(Z~=0);
    for i = randperm(N)
        m = m - (Z(i,:)~=0);
        chii = X(i,:) - Z(i,:)*V';

        for k = 1:K
            if m(k) == 0
                continue;
            end
            h(k) = h(k) - 2.^(N-i)*(Z(i,k)~=0);
            Vk = V(:,k);
            chii = chii + Z(i,k)*Vk';
            tinvsigma = sqrt(1./sigmaz^2+(Vk'*Vk)./sigmax^2);
            tmu = (chii*Vk)./(sigmax*tinvsigma)^2;
            pn0bp0 = sum(h==h(k))/(sum(h==h(k)+2.^(N-i))+1)*m(k)/(N-m(k))./(tinvsigma*sigmaz).*exp(0.5*(tmu*tinvsigma)^2); %sum(h==h(k))/(sum(h==h(k)+2.^(N-i))+1)*
            if algtype == 2 || algtype == 4            % MAP inference
                Z(i,k) = pn0bp0 > 1; % decide whether to augment the parameter space with one new latent variable (the real valued one)
                if Z(i,k) ~= 0       % and if do augment, pick its most probable value
                    Z(i,k) = tmu;
                    assert(tmu~=0);
                end
            else                                       % Gibbs sampling
                Z(i,k) = rand*(1+pn0bp0) > 1;
                if Z(i,k) ~= 0
                    Z(i,k) = randn./tinvsigma + tmu;
                end
            end
            h(k) = h(k) + 2.^(N-i)*(Z(i,k)~=0);
            chii = chii - Z(i,k)*Vk';
        end

        sglidx = find(m==0);
        chii = chii + Z(i,sglidx)*V(:,sglidx)';
        Z(i, sglidx) = 0;
        m = m + (Z(i,:)~=0);
        h(sglidx) = 0;

        rmXi = chii * chii';
        if algtype == 1 || algtype == 2                % TYPE-I
            kp = 0;
        else                                           % TYPE-II
            kp = min(numel(sglidx), poisstrunc);
        end
        loglhdk = zeros(1, poisstrunc);
        M = zeros(kp, kp);
        xp = zeros(1, kp);
        for newk = 1:kp
            tV = V(:,sglidx(1:newk));
            xp = sigmaz*(chii*tV);
            M = sigmaz^2*(tV'*tV);
            M(1:newk+1:newk*newk) = M(1:newk+1:newk*newk) + sigmax^2;
            loglhdk(newk) = -(D-newk)*log(sigmax) - 0.5*log(det(M)) + ...
                -0.5/sigmax^2*(rmXi - xp*(M\xp'));
        end
        for newk = kp+1:poisstrunc
            kh = newk-kp;
            loglhdk(newk) = log(integral(@fknuZnu, 0, Inf)); %, 'RelTol', 0.5e-10));
            loglhdk(newk) = loglhdk(newk) - kh*log(sqrt(2)*sigmav*sigmaz) - gammaln(kh/2);
        end
        logpk = [-D*log(sigmax)-0.5*rmXi/sigmax^2, ...
            cumsum(log(alphav/N)-log(1:poisstrunc)) + loglhdk];
        % check: log(sum(exp(logpk))) - (0.5/sigmax^2)*norm(X(setdiff(1:N,i),:)-Z(setdiff(1:N,i),:)*V', 'fro')^2
        if algtype == 2 || algtype == 4 % decide the number of new latent variables to be introduced
            [~, newk] = max(logpk);
            newk = newk - 1;
        else
            newk = logmnrnd(logpk(:))-1;
        end

        kh = max(0, newk-kp);
        newks(kh+1) = newks(kh+1) + 1;
        if kh == 0
            zetanew = 0;
            znewh = zeros(1,0);
        else
            mapzeta = argmaxfknuZnu();
            if algtype == 2 || algtype == 4
                zetanew = mapzeta;
                assert(zetanew>0);
                znewh = sqrt(zetanew/kh)./sigmav*ones(1,kh);
            else
                zetanew = slicesample(mapzeta, 1, 'logpdf', @(x)log(fknuZnu(x)));
                assert(zetanew>0);
                znewh = randn(1,kh);
                znewh = znewh.*(sqrt(zetanew)./sqrt(znewh*znewh')./sigmav);
            end
        end
        Z(i,K+1:K+kh) = znewh;
        m(K+1:K+kh) = 1;
        h(K+1:K+kh) = 2.^(N-i);

        kk = min(newk, kp);
        if kk == 0
            kpidx = [];
            znewp = zeros(1,0);
        else
            kpidx = sglidx(1:kk);
            tV = V(:,kpidx);
            invSigmap = (tV'*tV)./(sigmax^2+zetanew);
            invSigmap(1:kk+1:kk*kk) = invSigmap(1:kk+1:kk*kk) + 1/sigmaz^2;
            R = choll(invSigmap);
            mup = (R\(R'\(tV'*chii')))./(sigmax^2+zetanew);
            if algtype == 2 || algtype == 4
                znewp = mup';
            else
                znewp = (mup + R\randn(kk, 1))';
            end
            chii = chii - znewp*tV';
        end
        Z(i,kpidx) = znewp;
        m(kpidx) = 1;
        h(kpidx) = 2.^(N-i);

%         visualibp(Z,16,[0.2,0.8],false); pause(0.1);

        if kh > 0 % newk > kp
            tmu = kron(chii', (sigmav^2./(sigmax^2+zetanew))*znewh);
            if algtype == 2 || algtype == 4
                V(:,K+1:K+kh) = tmu;
            else
                tSigma = -(sigmav^4./(sigmax^2+zetanew))*(znewh'*znewh);
                tSigma(1:kh+1:kh*kh) = tSigma(1:kh+1:kh*kh) + sigmav^2;
                V(:,K+1:K+kh) = mvnrnd(tmu, tSigma);
            end
        end
%         K = K + kh;
        
        zcidx = ~any(Z);
        Z(:,zcidx) = [];
        V(:,zcidx) = [];
        m(zcidx) = [];
        h(zcidx) = [];
        K = size(Z,2);
    end
    assert(~any(~any(Z)));
%     newK = sum(newks(2:end).*(1:poisstrunc));
%     zcidx = ~any(Z(:,1:K-newK));
%     Z(:,zcidx) = [];
%     V(:,zcidx) = [];
%     K = size(Z,2);
    tElapsed = toc(tStart);
    if loopi > 0
        tfval = fobj(X, Z, V, alphav, sigmax, sigmaz, sigmav);
        fval = [fval; tfval];
        nfeapu = sum(Z~=0,2);
        maxnewkid = find(newks, 1, 'last');
        fprintf('[%d]: %d(%d), %s, [%s][%s], %.4f (%.2fs) | Z\n', loopi, ...
            sum(any(Z)), maxnewkid-1, num2str([mean(nfeapu),std(nfeapu(:),1)], ['%.2f',char(177),'%.2f']), ...
            num2str(newks(2:maxnewkid)./(sum(newks)-newks(1)), '%.2f '), ...
            num2str(newks(1:maxnewkid)./(sum(newks)-cumsum([0,newks(1:maxnewkid-1)])), '%.2f '), ...
            tfval, tElapsed);
    end
    etime = etime + tElapsed;
    
    % sampling V
    tStart = tic;
    tmp = Z./sigmax;
    invBj = tmp'*tmp;
    invBj(1:K+1:K*K) = invBj(1:K+1:K*K) + 1/sigmav^2;
    R = choll(invBj);
    b = (R\(R'\(Z'*X)))./sigmax^2;
    if algtype == 2 || algtype == 4
        V = b';
    else
        V = (b + R\randn(K, D))';
    end
    tElapsed = toc(tStart);
    if loopi > 0
        tfval = fobj(X, Z, V, alphav, sigmax, sigmaz, sigmav);
        fval = [fval; tfval];
        fprintf('[%d]: %s (%.2fs) | V\n', loopi, ... 
            num2str([mean(V(:)),std(V(:),1),tfval], ['%.2f',char(177),'%.2f',', %.4f']), tElapsed);
    end
    etime = etime + tElapsed;
end

    function p = fknuZnu(zeta)
        if kp == 0
            p = zeta.^(kh/2-1).*(zeta+sigmax^2).^(-D/2).* ...
                exp(-0.5*(zeta./(sigmav^2*sigmaz^2)+rmXi./(zeta+sigmax^2)));
        elseif kp == 1
            p = zeta.^(kh/2-1).*(zeta+sigmax^2).^(-(D-1)/2).*(M+zeta).^-0.5.* ...
                exp(-0.5*(zeta./(sigmav^2*sigmaz^2)+(rmXi-xp^2./(M+zeta))./(zeta+sigmax^2)));
        else
            p = zeta.^(kh/2-1).*(zeta+sigmax^2).^(-(D-kp)/2).*arrayfun(@fdettM,zeta).^-0.5.* ...
                exp(-0.5*(zeta./(sigmav^2*sigmaz^2)+(rmXi-arrayfun(@fxptMxp,zeta))./(zeta+sigmax^2)));
        end
        p(zeta<=0) = 0;
    end
    function optzeta = argmaxfknuZnu()
        diffknuZnu = [1, ((D+2-newk)*sigmav^2*sigmaz^2+2*sigmax^2), ...
            (sigmax^4+sigmav^2*sigmaz^2*((D+4-newk-kh)*sigmax^2-rmXi)), ...
            (2-kh)*sigmav^2*sigmaz^2*sigmax^4];
        zetas = roots(diffknuZnu);
        zetas = [1e-10; zetas(arrayfun(@isreal,zetas) & zetas >= 0)];
        [~, mapzetaid] = max(arrayfun(@fknuZnu, zetas));
        if kp > 0
            optzeta = fminsearch(@(x)-fknuZnu(x), zetas(mapzetaid), ...
                optimset('Display', 'off'));
        else
            optzeta = zetas(mapzetaid);
        end
    end
    function dettM = fdettM(zeta)
        if kp == 0
            dettM = 1;
            return;
        end
        if kp == 1
            dettM = M+zeta;
            return;
        end
        tM = M;
        tM(1:kp+1:kp*kp) = tM(1:kp+1:kp*kp) + zeta;
        dettM = det(tM);
    end
    function xptMxp = fxptMxp(zeta)
        if kp == 0
            xptMxp = 0;
            return;
        end
        if kp == 1
            xptMxp = xp^2./(M+zeta);
            return;
        end
        tM = M;
        tM(1:kp+1:kp*kp) = tM(1:kp+1:kp*kp) + zeta;
        xptMxp = xp*(tM\xp');
    end

end